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Abstract 

Previous models of disease spread involving delay have used basic SIR (susceptible - infectious - recovery) formulae and 
approaches. This paper demonstrates how time-varying SEIRS (S - exposed - I - R - S) models can be extended with delay 
to produce wave propagations that simulate periodic wave fronts of disease spread in the context of population 
movements. The model also takes into account the natural mortality associated with the disease spread. Understanding the 
delay of an infectious disease is critical when attempting to predict where and how fast the disease will propagate. We use 
cellular automata to model the delay and its effect on the spread of infectious diseases where population movement occurs. 
We illustrate an approach using wavelet transform analysis to understand the impact of the delay on the spread of 
infectious diseases. The results indicate that including delay provides novel ways to understand the effects of migration and 
population movement on disease spread. 
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Introduction 

There is very little understanding of how disease spread is 
affected by population migration. When modeling the propagation 
of disease, it is important to model not only individual infections 
but also the spread of infection from populations of infected 
individuals to other populations [1]. In particular, if the disease is 
fatal, what matters is the likelihood of the disease spreading to new 
individuals before the currently infected individuals die. Mathe- 
matical models have been developed for use in the dynamic 
modeling of disease spread [2,3], such as the SIR model, which 
represents the state changes among the members of the susceptible 
(6), infected (7) and recovered [R) populations. The goal of such 
models is to understand the disease spread patterns and to predict 
the outcomes of introducing public health intervention measures 
to minimize the spread of disease. Such models must take into 
account the difference in time between becoming exposed {E) to 
infection and becoming infected (hence, SEIR) as well as between 
becoming infected and dying ('delays'). A wide variety of such 
approaches of modeling disease with delay have been formulated. 
Some of these delay concepts are strongly parametric (specific 
distributions), while others are nonparametric (general distribu- 
tions). However, there is very little understanding of how the SEIR 
models are affected by delay [4]. 

A delay in many population models considers that the 
transmission dynamic behavior of the disease at time t can 
destabilize the equilibrium [5]. There are two types of delay: a 
discrete or fixed delay and continuous or distributed delay [6] . In 
models with a discrete or fixed delay, the dynamic behavior of the 
disease depends on the state at time t — T, where t is a fixed 



constant [7]. For example, the number of newborns at time t 
depends on the state of population at time t — T, where t the period 
of pregnancy. In the case of a continuous or distributed delay, the 
dynamic behavior of the model at time t depends on the states 
during the entire period prior to time t. Both types of delays 
produce outputs at the end of each time step that can be 
interpreted as a signal, or function, of the model. This signal or 
function can be analyzed further to extract the features, such as 
frequency and amplitude, which can be mapped onto interpret- 
able properties of the model. Starting a model with two different 
sets of parameter values can lead to two different signals over time, 
which, upon further analysis, can be mapped back to the effects 
that the parameter values have on the model. Such signals can 
include oscillations and waves if periodicity (a function that repeats 
its values at regular intervals) emerges from the model. The 
periodicity can itself be analyzed using wavelet or Fourier 
transforms. 

The advantage of the wavelet transform of a signal over the 
Fourier transform is the ability of wavelet transforms to capture 
local information (i.e., feature changes in space or time), whereas 
Fourier transforms capture global features (e.g., harmonic 
components of the entire signal). That is, wavelet transforms are 
localized in both time and space, and can be useful for identifying 
parts of the signal that are local in terms of the geographical 
location. Wavelet transforms have many uses in the areas of wave 
propagation, data compression, signal processing, image process- 
ing and pattern recognition [8,9]. In this paper, we demonstrate 
how cellular automata supplemented with wavelet transforms can 
be used to implement a SEIRS model that attempts to simulate the 
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spread of disease locally in both time and space, where space refers 
to the delay from exposed individuals to infectious individuals. 

Description of the Model 

In this section, we introduce cellular automata (CA) for 
modeling geographical distributions. CA has a significant role in 
epidemic modeling because each individual or cell, or a small 
region of space, updates itself independently, allowing for the 
concurrent development of several epidemic spatial clusters. The 
new state of a cell is based on its current state and the current 
states of its surrounding cells as well as on some shared rules of 
change. 

The SEIRS model used here is specified as follows. Let Sij be 
the number of susceptible individuals / in the population at some 
patch j at time t. Let Eij represent the number of susceptible 
individuals who become exposed (infected but not infectious) [10] 
at that patch. An individual remains in the exposed class or the 
exposed state for a certain latent period before becoming infected 
[11]. Let lij denote the number infected at that patch and Rfj 
denote the number recovered at that patch. If > 0 is the 
effective contact rate per individual per unit of time at a patch, we 
can use the Law of Mass Action assumption that the rate of change 
from being susceptible to being exposed and finally to natural 
death is proportional to the population number in each state, with 
a rate constant jd. The assumption is that t is the delay time (days), 
which is a constant. Adopting the approach of Yan and Liu [12], 
the probability that an individual survives the latent period to 
become infected at [t — T,t] is e~^K Additionally, the number of 
susceptible individuals who becomes exposed at time t — T is 
PijSij{t — z)Iij{t — z)/Nij{t — %). The difference between 
PijSij{t)Iij{t) and PijSij{t — x)Iij{t — x) is that former is the exposed 
state and the latter is the infected state, which leads to a 

deterministic model of d a ^^^^-^r-)-^. — —e~^^. 

Nij(t-T) 

The task for the wavelet transform is to identify, from the 
composite signal formed from the outputs of the model at each 
time step (to be described in more detail below), any local temporal 
and spatial periodicities that emerge. 

Rules for Disease Spread 

The rules described below [13] determine the state transitions of 
individual cells in the CA for the SEIRS model at each time step 
and will incorporate other probabilistic parameters. 

• A cell changes its state from susceptible to exposed (S^E) when it 
comes in contact with an infected cell within its defined 
neighborhood. 




Figure 1. Movement with all exposed class in SEIRS Model. 

Depicts the flow chart of the SEIRS model. 
doi:10.1371/joumal.pone.0098288.g001 



• A cell changes its state from exposed to infectious (E^I) after 
being in the state E for a given e, which is the transition time. 

• The state of the cell changes from infectious to recovery (I^R) 
after being in the state / for a given time a. In our model, we 
assume that for every cell, the amount of time a elapsed from 
the infectious cell state to the recovery cell state. In state /, the cells 
are capable of passing on the infection to neighboring cells. 

• The cell remains in state (R^S), signifying complete recovery 
for some time 3, which is the transition time between the 
recovery state and the susceptible state, as shown in Fig. 1. 

These CA rules can be described by a concrete expression using 
the classic SEIRS models based on differential equations. The set 
of ordinary differential equations corresponding to the CA model 
is described in the following. 

dlii ^ , Sij{t — T)Iij{t — T) , _ 

N,(t-^ -(M+o<)/,(0 + .£,(0 (3) 

dR 

-^=OClyit)-itl + d)Rijit) (4) 

where k = e^'" 

Nij(t) = Sijit) + Eyit) + lyit) + Rijit) 

It is assumed that the population sizes of each patch are 
identical and remain constant. This corresponds to assuming that 
the population birth and death rates, denoted by b and fi 
respectively. Individuals are born susceptible and can acquire 
infection from infective individuals, at which point they enter the 
exposed state. After a latent period, lasting an average of 1/e time 
units, individuals are infectious for an average of 1/a time units. 
Since the population is subject to a disease-independent mortality 
rate, /i, the mortality-corrected average duration of infectiousness 
is !/(/! + a), and the exposed individuals who will become 
infectives is + This formulation of the model assumes that 
there is migration rate, m into exposed class to ascertain the period 
at which all exposed class migrate to infected class. 

As mentioned above, the infection is modeled as the contact 
between an infected individual and susceptible individuals in the 
neighborhood, either as a deterministic or as a probabilistic 
process, depending on how the conditions on the arcs of the model 
in Figure 1 are realized. For our experiments below, we use 
probabilistic conditions. 

Simulation Scenario 

In this paper, r is the radius used for the neighborhood around 
the cardinality point within which the nearest cell is infected. The 
order Moore neighborhood was defined as the 8 nearest 
neighbors, i.e., r=l. The degree of infectiousness was used to 
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Table 1. Simulation Protocol. 



Events 


Model Prediction 


CA model 


Environment 


Spatial (landscape) 


Spatial (landscape) 


Susceptible 


9980 people 


9980 people 


Exposed 


Nil 


From model prediction 


Infected 


20 people 


20 people 


Recovery 


Nil 


From model prediction 


Natural death rate 


0.0005 yr^ 


0.0005 yr^ 


Incubation period, t 


7 days 


7 days 


No. of simulation 


365 iterations 


365 iterations 


Infected Factor, p 


0.025 


0.25 



doi:1 0.1 371/journal.pone.0098288.t001 



assess the probabilistic distribution that the current cell is infected 
in the next time step. The disease will propagate through the 
landscape based on a set of probabilities of state transitions [14]. 
Before an infected cell disperses (migrates), it updates indepen- 
dently, based on the variable location of the neighboring cells, with 
probability p 



2 * (infected parameter ) 
" (2r+ 1)2-1 ' 



(5) 



where (2r+ 1) — 1 is the cardinality of the neighboring cell. 

Each infected individual can enter the cell from the cardinality 
of the neighboring cell with the same probability p, whereas there 
are a number of new infected that arrive in the cell from the 
neighboring cells. Finally, the state of the cell after a complete 
transition at time t is updated by the iteration number (time), that 
is, t is increased by 1 to become time (/+ 1). 

Stability Analysis 

An important equilibrium point for any disease model is the 
disease-free equilibrium (DFE). The stability of the DFE is 
especially important because it determines whether a disease is 
capable of attacking an entire population beyond the prior 
expectation (epidemic or pandemic) [15]. The reproduction 
number Rq is a threshold value or number that determines the 



stability of the DFE [16]. The reproduction number is the 
expected number of secondary cases produced by a typical 
infection in a completely susceptible population [17]. If i?o > 1? an 
epidemic occurs; if Ro<l, an epidemic does not occur; and if 
Ro = \, a change of stability occurs. Rq can be calculated as 
Rq = p{FV~^\ the spectral radius (p) of the next generation 
matrix FV~^ . The {ij) is the number of entries of matrix FV~^ , 
which represents the number of new infections in compartment / 
due to an infected individual being introduced into compartment j 
[17,18]. 

Let Fi{x) be the rate at which newly infected individuals 
(transmitted rate) enter compartment and let Vi = Vf — V^^ 
denote the transfer of individuals out of ( Vj~ ) and into ( V^^ ) the 
i^^ compartment. With this interpretation, we write a matrix F 
that defines the rate of new infections in different compartments, 
differentiated with respect to E and /, and then evaluated at the 
disease-free equilibrium. An equilibrium solution with 




Figure 2. A dynamic spread process. Depict the geographical 
distribution of the propagation of disease spread clusters for 40. 
doi:1 0.1 371 /journal.pone.0098288.g002 



Figure 3. A dynamic spread process. Depict the geographical 
distribution of the propagation of disease spread clusters for 65 
simulations. 

doi:10.1371/journal.pone.0098288.g003 



PLOS ONE I www.plosone.org 



3 



June 2014 | Volume 9 | Issue 6 | e98288 



Wavelet Propagation Distribution of SEIRS Modeling with Delay 





E = I = R = 0 and S = So, where *So is any positive solution of 
b = iiSo. Therefore, this will be locally stable and hence a disease 
free equilibrium, So = b/ ji. We assume that, this to be the case, 
evaluating the derivatives of F and V a^t S = Sq, E = 1 = R = 0, we 
get the following expressions for F, V and the spectral radius 
respectively. 



0 f^yk 



Now, we also write V that defines the rate of transfer of infected 
from one compartment to another (F/= Vf — F^-^). 

/;U + e 0 \ 
^= ^ (7) 



FV-' = 



(;U + £)(^ + a) 



(8) 



Figure 4. A dynamic spread process without exposed state. 

Depicts how the propagation of the disease spread clusters proceeds in 
geographical distribution using 365 simulations. The model indicates 
that there is increasing volatility in the susceptible population after 240 
days due to outgoing moves from the susceptible to the exposed 
population (Figure 5A). 
doi:10.1371/journal.pone.0098288.g004 



The reproduction number Rq is given by the dominant eigenvalue 
(or spectral radius) of FV~^ . 
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Figure 5. Snapshot of the susceptible population. The model indicates that there is increasing volatility in the susceptible population after 240 
days due to outgoing moves from the susceptible to the exposed population (Figure 5A). The spread of the disease in the discrete transform 
exhibited peaks at 230 days and 360 days (Figure 5B), as indicated by the arrows. Slower and faster frequencies were found at the initial stages as well 
as the final stages at 10 days and 360 days, as indicated by the arrows. The blue regions denote the probability densities. For example, at 330 days, 
there was an indication of a slower frequency that resulted in a faster spread of the disease (Figure 5C) due to the large coefficients in the continuous 
wavelet transform. 

doi:1 0.1 371 /journal. pone.0098288.g005 
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Figure 6. Snapshot of the exposed population. The model indicates that there is an increased volatility in the exposed population after 200 
days, due to the outgoing moves from the exposed population to the infectious population (Figure 6A), and an increased volatility in the infectious 
population after 240 days, due to the outgoing moves from the infectious population to the recovery population (Figure 7A). The spread of the 
disease in the discrete transform reached a peak at 70 days due to the large coefficients (Figure 6B). In the continuous wavelet transform at 72 days, 
there were no exposed individuals (Figure 6C) because they moved to the infectious stage. 
doi:10.1371/journal.pone.0098288.g006 
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Figure 7. Snapshot of the exposed population. The spread of the disease in the discrete transform peaked at 240 days (Figure 7B), and at 240 
days and 280 days (Figure 7C) due to large coefficients. The range of frequencies used in averaging is indicated by the arrow at 240 days, which 
corresponds to the peak of the disease spread (Figure 7C). 
doi:1 0.1 371 /journal. pone.0098288.g007 
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Figure 8. Snapshot of the infected population. The model indicates that there is increasing volatility in the recovery population after 240 days 
due to the outgoing moves from the recovery population to the natural death population (Figure 8A). The spread of the disease in the discrete 
transform peaked at 240 days (Figure 8B) and at 230, 250 and 340 days due to large coefficients (Figure 8C). The range of frequencies used in 
averaging is indicated by the arrow at 250 and 362 days, which correspond to the peaks of the disease spread (Figure 8C). 
doi:1 0.1 371 /journal. pone.0098288.g008 
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Figure 9. Snapshot of the recovery population. An increase in natural deaths was observed after 50 days due to the outgoing moves from the 
infectious and recovery populations to the natural death population (Figure 9A). The range of frequencies used in averaging is indicated by the arrow 
at 362 days, which corresponds to the peak of the disease spread (Figure 9C). 
doi:1 0.1 371 /journal. pone.0098288.g009 
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Figure 10. Sinusoidal response of the exposed population in 3 
dimensions. The CWT coefficients are large at scales near the 
frequencies of the exposed waves and clearly depict the sinusoidal 
pattern in the CWT coefficients at these scales for the exposed 
population. 

doi:10.1371/journal.pone.0098288.g010 

Then, 
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CWT of Exposed. 

35 r : : : ; 




5- 



qI \ \ i \ i \ \ i 

0 10 20 30 40 50 60 70 80 

Days 

Figure 1 1 . Sinusoidal response of the exposed population. Plots 
the same transform from Figure 10 but with different angle for better 
visualization, depicts the time period at which all the individuals who 
were exposed to disease moved to infectious population. 
doi:10.1371/journal.pone.0098288.g011 

In that case, we consider a wavelet function ij/, such that 



Ro- 



(n-\-a){H + e) 



(9) 



Ro = p{FV~^), where p denotes the spectral radius. 
If 7?o<l, then the DFE is globally asymptotically stable; if 
Rq>\, then the DFE is unstable [16]. 

Why Use Wavelet Analysis? 

The wavelet transform is a function that is an improved 
technique of implementing the Fourier transformation. In the 
time-frequency analysis of a set of data, the classical Fourier 
transform analysis is inadequate due to a lack of any local 
information contained in the data. The use of wavelet analysis 
describes the pattern, trends, and the structures that might be 
overlooked in the raw data. The usefulness of wavelets in data 
analysis is very clear, particularly in the field of statistics, where 
large and cumbersome data sets are prevalent. Wavelet analysis is 
used here to accomplish the following tasks: 

a) to analyze the frequency of the spread of the disease; and 

b) to understand the impact of the delay in the spread of an 
infectious disease. 

These two tasks will be accomplished by analyzing the 
frequency and interpreting the exposed as a function of the delay. 
Let a be certain class of functions and /o/i5/2v/« be the 

GO 

frequencies of a simple function such that each f{x)= J] cinfn{x) 

n = 0 

for some coefficients a„ . We will consider a function in the Laplace 
transform L2(R), i.e.. 



^(x)P<« 



(10) 



jeZ keZ 



(11) 



where dj^k are the wavelet coefficients and ^jj^{x) = i'^^{7Jx — k) 
are the translated and scaled version of the wavelet xj/. 

Let s represent each compartment and ^ be the wavelet. Then, 
the wavelet coefficient of s at scale a and position b is defined by: 



s{t)^\l/l—^\dt 
yja \ a 



(12) 



Because s{t) is discrete, we will use a piecewise constant 
interpolation of the s{k) values, where k= \ to length[s). 

Simulation Setup 

Cellular automata are found to have received extensive 
academic study as one of the mathematical tools for successfully 
modeling the spread of diseases. Because there are a number of 
cellular automata used in the Matlab programming environment 
[19], Matlab was chosen as the implementation tool here. 

Experiment and Results 

We constructed cellular automata in Matlab to implement the 
model described by equations (l)-(4) stated in the previous section. 
The experiments were performed on a 100x100 cell spaces. We 
assumed that twenty infected individuals were introduced into the 
population of 100x100 cells, which corresponds to ten-thousand 
individuals. Table 1 summarizes the parameter values used. Blue, 
light green, red and light blue squares correspond to the values of 
0, 0.5, 1 and 0.1, respectively, with the values of 0, 0.5, 1 and 0.1, 
denoting susceptible, exposed, and infected and recovery, respec- 
tively. A total of 365 iterations were simulated, with one iteration 
representing one day. The variables and parameters are shown in 
Table 1. 
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Figure 2 and 3 depict the geographical distribution of the 
propagation of disease spread clusters for 40 and 65 simulations, 
respectively. 

The colors dark blue, green, red and light blue represent healthy 
individuals (susceptible), infected individuals (exposed), infectious 
individuals and recovery individuals, respectively. We can see the 
breakdown of an initially homogenous spread-of- disease pattern. 
As the phase separation progresses, a persistent compact spread of 
disease is formed, as shown in Figure 2. Figure 3 shows the further 
progression of the phase separation, in which a persistent compact 
spread of disease is formed (infectious), surrounded by exposed 
individual populations. 

Figure 4 depicts how the propagation of the disease spread 
clusters proceeds in geographical distribution using 365 simula- 
tions. The model indicates that there is increasing volatility in the 
susceptible population after 240 days due to outgoing moves from 
the susceptible to the exposed population (Figure 5A). 

The colors of dark blue, red and light blue represent healthy 
individuals (susceptible), infectious individuals and recovery 
individuals, respectively. As the phase separation proceeds, a 
persistent compact spread of disease is formed (infectious and 
recovery) within the population. The spread of the disease in the 
discrete transform exhibited peaks at 230 days and 360 days 
(Figure 5B), as indicated by the arrows. Slower and faster 
frequencies were found at the initial stages as well as the final 
stages at 10 days and 360 days, as indicated by the arrows. The 
blue regions denote the probability densities. For example, at 330 
days, there was an indication of a slower frequency that resulted in 
a faster spread of the disease (Figure 5G) due to the large 
coefficients in the continuous wavelet transform. 

The model indicates that there is an increased volatility in the 
exposed population after 200 days, due to the outgoing moves 
from the exposed population to the infectious population 
(Figure 6A), and an increased volatility in the infectious population 
after 240 days, due to the outgoing moves from the infectious 
population to the recovery population (Figure 7 A). 

The spread of the disease in the discrete transform reached a 
peak at 70 days due to the large coefficients (Figure 6B). In the 
continuous wavelet transform at 72 days, there were no exposed 
individuals (Figure 6C) because they moved to the infectious stage. 
The spread of the disease in the discrete transform peaked at 240 
days (Figure 7B), and at 240 days and 280 days (Figure 7C) due to 
large coefficients. The range of frequencies used in averaging is 
indicated by the arrow at 240 days, which corresponds to the peak 
of the disease spread (Figure 7C). 

The model indicates that there is increasing volatility in 
recovery after 240 days due to the outgoing moves from the 
infectious population to the recovery population (Figure 8A). In 
addition, an increase in natural deaths was observed after 50 days 
due to the outgoing moves from the infectious and recovery 
populations to the natural death population (Figure 9A). 

The model indicates that there is increasing volatility in the 
recovery population after 240 days due to the outgoing moves 
from the recovery population to the natural death population 
(Figure 8A). The spread of the disease in the discrete transform 
peaked at 240 days (Figure 8B) and at 230, 250 and 340 days due 
to large coefficients (Figure 8C). The range of frequencies used in 
averaging is indicated by the arrow at 250 and 362 days, which 
correspond to the peaks of the disease spread (Figure 8C). The 
spread of the disease in the discrete transform were at a peak at 
240 days (Figure 9B). The range of frequencies used in averaging is 
indicated by the arrow at 362 days, which corresponds to the peak 
of the disease spread (Figure 9C). 



As shown in Figure 10, the CWT coefficients are large at scales 
near the frequencies of the exposed waves and clearly depict the 
sinusoidal pattern in the CWT coefficients at these scales for the 
exposed population. Figure 1 1 plots the same transform from a 
different angle for better visualization. 

The sinusoidal wave amplitude is the height of the crest, and the 
frequency is the number of oscillations per second. Hence, the 
amplitude remains the same for any change in frequency. The 
maximum and minimum are both at 225 Hz, which demonstrates 
that the spread of the disease in the exposed population was at the 
peak state (Figure 10). In this plot, the value of each (x, y) 
coordinate represents the strength of the spread of the disease 
between the coordinates. The strong interaction between the (x,y) 
coordinates is 0.026 Hz (Figure 1 1). As seen in Fig. 6C, the arrow 
at 70 days indicated the faster frequency, which resulted in a 
slower spread of the disease. 

Conclusion 

We have developed a mathematical disease propagation model 
of a susceptible-infected-recovery (SIR) type extended to a four 
compartmental epidemiological model with a susceptible-exposed- 
infectious-recovery-susceptible (SEIRS) with exposed class moved 
as the spread of an infectious disease in a population with a delay. 
The discrete and continuous transforms demonstrate where and 
how the propagation of disease was dominant in the population. 
The interactions between the infected and the susceptible lead to 
an exposed stage of the disease, which agrees with global dynamics 
behaviors for a new delay SEIR epidemic disease model proposed 
by Meng and Chen et al. [4]. Wavelets are very important for 
detecting abrupt changes in the distribution of disease spread [20] . 
These abrupt changes occur from one compartment to another 
(susceptible, exposed, infectious, recovery, and susceptible). These 
changes produced relatively large wavelet coefficients centered on 
the discontinuity at all scales. The location of the discontinuity 
based on the CWT coefficients is obtained at the smallest scales. 
This selection for shorter infection period was also reported in 
[13]. Note that, due to the stochastic nature of the model 
parameters, the introduction of the dead period is only a transient 
phenomenon. As shown in Figure 7, the number of individuals 
who were infectious to the disease was much higher at the 
beginning of the spread of the disease. An infection will only cause 
a new infection if the individual is still exposed at the time of 
infection. Alimadad and Dabbaghian et al. [21] also reported that 
an individual becomes exposed to an infection when a randomly 
chosen neighbor individual is susceptible to it. The dominant 
order of the frequencies indicates that a delay in the transition 
between exposed to infectious states occurs after the peak points in 
the different time intervals. The detection of these discontinuities 
(delays) was associated with the speed of the spread of the disease 
and is detectable in the frequency and phase of the CWT. The 
disease stopped spreading in the exposed state at 72 days 
(Figure 1 1) but continued to spread in the susceptible, infectious 
and recovery populations after 240 days. The spread of the 
recovery population continued to persist, as shown in Figure 4, 5, 
7, 8 and 9. This behavior implies that: a) the higher the frequency, 
the lower the spread of the disease, and b) the faster the spread of 
the disease, the slower the frequency. The discrete and continuous 
transforms demonstrate where and how the propagation of disease 
was dominant in the population. Therefore, we have demonstrat- 
ed how to model the delay through exposure and the effects of the 
delay on a movement of exposed state to infectious population. 
Finally, the proposed model will be tested against real world data 
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to test its fitness and to evaluate its accuracy as well as to determine 
the utility of the cellular automata model. 
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